#Anti-Muslim Policy Preferences and Boundaries
#of American Identity Across Partisanship
#Additional Analyses

#Clear R
rm(list=ls())

#Packages
library(tidyverse)
library(haven)
library(margins)
library(stargazer)
library(ggthemes)
library(broom)
library(car)
library(magrittr)
library(ggridges)
library(survey)
library(descr)
library(reshape2) 
library(anesrake)
library(rstatix)
library(ltm) 
library(data.table) 
library(misty)
library(sensemakr)
library(lavaan)



#################
# Load Datasets #
#################

#2018 CCES
CCES<-read_dta("2018_10_CCES_final_weight.dta")
glimpse(CCES)

#2019 6 Lucid
Lucid_2019<-read_dta("2019_6_Lucid_final_weight.dta")
glimpse(Lucid_2019)

#2020 8 Lucid
Lucid_2020<-read_dta("2020_8_Lucid_final_weight.dta")
glimpse(Lucid_2020)


#################################################
# Examining potential acquiescence bias for key IV #
#################################################

#Filtering out respondents who selected all 1s and all 5
CCES_AQ1<-CCES %>% filter(Amid>0) #14 respondents or 1.9%
CCES_AQ2<-CCES %>% filter(Amid<1) #187 respondents or 25.9%
CCES_AQ<-CCES %>% filter(Amid>0 & Amid<1)


Lucid_2019_AQ1<-Lucid_2019 %>% filter(Amid>0) #34 or 0.9%
Lucid_2019_AQ2<-Lucid_2019 %>% filter(Amid<1) #844 or 22.6%
Lucid_2019_AQ3<-Lucid_2019 %>% filter(Amid>0 & Amid<1)

Lucid_2020_AQ1<-Lucid_2020 %>% filter(Amid>0) #27 or 1.4%
Lucid_2020_AQ2<-Lucid_2020 %>% filter(Amid<1) #590 or 31.4%
Lucid_2020_AQ<-Lucid_2020 %>% filter(Amid>0 & Amid<1)



############
# Modeling #
# CCES     #
############

#(i) How strong would a particular confounder (or group of confounders) have to be to change our conclusions? 
#(ii) In a worst case scenario, how vulnerable is our result to many or all unobserved confounders acting together, possibly non-linearly? 
#(iii) Are the confounders that would alter our conclusions plausible, or at least how strong would they have to be relative to observed covariates?


#Muslim Ban
summary(model_banmuslims_cces<-lm(banmuslims ~Amid + female + white + income + educ+  agegroup +
                                    attention + ideology_7pt + independent_other + democrat + 
                                    presapp , weight=teamweight, data= CCES))


banmuslims_cces.sensitivity<-sensemakr(model = model_banmuslims_cces, 
                                       treatment = "Amid",
                                       benchmark_covariates = "presapp",
                                       kd = 1:3,
                                       ky = 1:3, 
                                       q = 1,
                                       alpha = 0.05, 
                                       reduce = TRUE)


ovb_minimal_reporting(banmuslims_cces.sensitivity, format = "latex")


banmuslims_cces.sensitivity


plot(banmuslims_cces.sensitivity)




#Registry 
summary(model_registry_cces<-lm(registry ~Amid + female + white + income + educ+  agegroup +
                                  attention + ideology_7pt + independent_other + democrat + 
                                  presapp , weight=teamweight, data= CCES))



registry_cces.sensitivity<-sensemakr(model = model_registry_cces, 
                                     treatment = "Amid",
                                     benchmark_covariates = "presapp",
                                     kd = 1:3,
                                     ky = 1:3, 
                                     q = 1,
                                     alpha = 0.05, 
                                     reduce = TRUE)

registry_cces.sensitivity


ovb_minimal_reporting(registry_cces.sensitivity, format = "latex")




#Sharia 
summary(model_sharia_cces<-lm(sharia ~Amid + female + white + income + educ+  agegroup +
                                attention + ideology_7pt + independent_other + democrat + 
                                presapp , weight=teamweight, data= CCES))



sharia_cces.sensitivity<-sensemakr(model = model_sharia_cces, 
                                   treatment = "Amid",
                                   benchmark_covariates = "presapp",
                                   kd = 1:3,
                                   ky = 1:3, 
                                   q = 1,
                                   alpha = 0.05, 
                                   reduce = TRUE)

sharia_cces.sensitivity


ovb_minimal_reporting(sharia_cces.sensitivity, format = "latex")



##############
# Modeling   #
# 2019 Lucid #
##############

#Muslim Ban
summary(model_banmuslims_2019lucid<-lm(banmuslims ~Amid + female + white + income + educ+  agegroup +
                                         attention + ideology_7pt + independent_other + democrat + 
                                         presapp , weight=weight, data= Lucid_2019))



banmuslims_2019lucid.sensitivity<-sensemakr(model = model_banmuslims_2019lucid, 
                                            treatment = "Amid",
                                            benchmark_covariates = "presapp",
                                            kd = 1:3,
                                            ky = 1:3, 
                                            q = 1,
                                            alpha = 0.05, 
                                            reduce = TRUE)


ovb_minimal_reporting(banmuslims_2019lucid.sensitivity, format = "latex")



#Patrol
summary(model_patrol_2019lucid<-lm(patrol ~Amid + female + white + income + educ+  agegroup +
                                     attention + ideology_7pt + independent_other + democrat + 
                                     presapp , weight=weight, data= Lucid_2019))




patrol_2019lucid.sensitivity<-sensemakr(model = model_patrol_2019lucid, 
                                        treatment = "Amid",
                                        benchmark_covariates = "presapp",
                                        kd = 1:3,
                                        ky = 1:3, 
                                        q = 1,
                                        alpha = 0.05, 
                                        reduce = TRUE)


ovb_minimal_reporting(patrol_2019lucid.sensitivity, format = "latex")



#Sharia
summary(model_sharia_2019lucid<-lm(sharia ~Amid + female + white + income + educ+  agegroup +
                                     attention + ideology_7pt + independent_other + democrat + 
                                     presapp , weight=weight, data= Lucid_2019))




sharia_2019lucid.sensitivity<-sensemakr(model = model_sharia_2019lucid, 
                                        treatment = "Amid",
                                        benchmark_covariates = "presapp",
                                        kd = 1:3,
                                        ky = 1:3, 
                                        q = 1,
                                        alpha = 0.05, 
                                        reduce = TRUE)


ovb_minimal_reporting(sharia_2019lucid.sensitivity, format = "latex")



#restrictmosques
summary(model_restrictmosques_2019lucid<-lm(restrictmosques ~Amid + female + white + income + educ+  agegroup +
                                              attention + ideology_7pt + independent_other + democrat + 
                                              presapp , weight=weight, data= Lucid_2019))




restrictmosques_2019lucid.sensitivity<-sensemakr(model = model_restrictmosques_2019lucid, 
                                                 treatment = "Amid",
                                                 benchmark_covariates = "presapp",
                                                 kd = 1:3,
                                                 ky = 1:3, 
                                                 q = 1,
                                                 alpha = 0.05, 
                                                 reduce = TRUE)


ovb_minimal_reporting(restrictmosques_2019lucid.sensitivity, format = "latex")



##############
# Modeling   #
# 2020 Lucid #
##############


#Muslim Ban
summary(model_banmuslims_2020lucid<-lm(banmuslims ~Amid + female + white + income + educ+  agegroup +
                                         attention + ideology_7pt + independent_other + democrat + 
                                         presapp , weight=weight, data= Lucid_2020))




banmuslims_2020lucid.sensitivity<-sensemakr(model = model_banmuslims_2020lucid, 
                                            treatment = "Amid",
                                            benchmark_covariates = "presapp",
                                            kd = 1:3,
                                            ky = 1:3, 
                                            q = 1,
                                            alpha = 0.05, 
                                            reduce = TRUE)


ovb_minimal_reporting(banmuslims_2020lucid.sensitivity, format = "latex")





#noweapon
summary(model_noweapon_2020lucid<-lm(noweapon ~Amid + female + white + income + educ+  agegroup +
                                       attention + ideology_7pt + independent_other + democrat + 
                                       presapp , weight=weight, data= Lucid_2020))




noweapon_2020lucid.sensitivity<-sensemakr(model = model_noweapon_2020lucid, 
                                          treatment = "Amid",
                                          benchmark_covariates = "presapp",
                                          kd = 1:3,
                                          ky = 1:3, 
                                          q = 1,
                                          alpha = 0.05, 
                                          reduce = TRUE)


ovb_minimal_reporting(noweapon_2020lucid.sensitivity, format = "latex")




#sharia
summary(model_sharia_2020lucid<-lm(sharia ~Amid + female + white + income + educ+  agegroup +
                                     attention + ideology_7pt + independent_other + democrat + 
                                     presapp , weight=weight, data= Lucid_2020))




sharia_2020lucid.sensitivity<-sensemakr(model = model_sharia_2020lucid, 
                                        treatment = "Amid",
                                        benchmark_covariates = "presapp",
                                        kd = 1:3,
                                        ky = 1:3, 
                                        q = 1,
                                        alpha = 0.05, 
                                        reduce = TRUE)


ovb_minimal_reporting(sharia_2020lucid.sensitivity, format = "latex")





#restrictmosques
summary(model_restrictmosques_2020lucid<-lm(restrictmosques ~Amid + female + white + income + educ+  agegroup +
                                              attention + ideology_7pt + independent_other + democrat + 
                                              presapp , weight=weight, data= Lucid_2020))




restrictmosques_2020lucid.sensitivity<-sensemakr(model = model_restrictmosques_2020lucid, 
                                                 treatment = "Amid",
                                                 benchmark_covariates = "presapp",
                                                 kd = 1:3,
                                                 ky = 1:3, 
                                                 q = 1,
                                                 alpha = 0.05, 
                                                 reduce = TRUE)


ovb_minimal_reporting(restrictmosques_2020lucid.sensitivity, format = "latex")




#########################
#Histogram of AMID by PID
#########################

#CCES
CCES<-CCES %>% mutate(pid=case_when(democrat==1~0, independent_other==1~1, republican==1~2))

# Take just American Identity Variable
iv_dat <- CCES[, c("Amid", "pid")]
iv_dat <- na.omit(iv_dat) # Omit missing data 

facetedbars <- ggplot(data = iv_dat, aes(Amid))  +
  geom_histogram(aes(y =..density..*.1), binwidth =.1) +
  scale_y_continuous("Percent", labels = scales::percent) +
  scale_x_continuous("American Identity (low to high)", limits=c(-0.05,1.14), breaks=seq( 0, 1, .1)) +
  scale_colour_identity() + 
  facet_wrap( ~ factor(pid, labels=c("Democrat", "Independent", "Republican")), 
              nrow=1) 
facetedbars + aes(fill = as.factor(-pid)) + theme_minimal() +theme(legend.position = "none") 
ggsave("cces_2018_amid_pid.png", 
       width = 12, height = 6)



#Lucid 2019
Lucid_2019<-Lucid_2019 %>% mutate(pid=case_when(democrat==1~0, 
                                                            independent_other==1~1,
                                                            republican==1~2))

# Take just American Identity Variable
iv_dat <- Lucid_2019[, c("Amid", "pid")]
iv_dat <- na.omit(iv_dat) # Omit missing data 

facetedbars <- ggplot(data = iv_dat, aes(Amid))  +
  geom_histogram(aes(y =..density..*.1), binwidth =.1) +
  scale_y_continuous("Percent", labels = scales::percent) +
  scale_x_continuous("American Identity (low to high)", limits=c(-0.05,1.14), breaks=seq( 0, 1, .1)) +
  scale_colour_identity() + 
  facet_wrap( ~ factor(pid, labels=c("Democrat", "Independent", "Republican")), 
              nrow=1) 
facetedbars + aes(fill = as.factor(-pid)) + theme_minimal() +theme(legend.position = "none") 
ggsave("lucid_2019_amid_pid.png", 
       width = 12, height = 6)





#Lucid 2020
Lucid_2020<-Lucid_2020 %>% mutate(pid=case_when(democrat==1~0, 
                                                            independent_other==1~1,
                                                            republican==1~2))

# Take just American Identity Variable
iv_dat <- Lucid_2020[, c("Amid", "pid")]
iv_dat <- na.omit(iv_dat) # Omit missing data 

facetedbars <- ggplot(data = iv_dat, aes(Amid))  +
  geom_histogram(aes(y =..density..*.1), binwidth =.1) +
  scale_y_continuous("Percent", labels = scales::percent) +
  scale_x_continuous("American Identity (low to high)", limits=c(-0.05,1.14), breaks=seq( 0, 1, .1)) +
  scale_colour_identity() + 
  facet_wrap( ~ factor(pid, labels=c("Democrat", "Independent", "Republican")), 
              nrow=1) 
facetedbars + aes(fill = as.factor(-pid)) + theme_minimal() +theme(legend.position = "none") 
ggsave("lucid_2020_amid_pid.png", 
       width = 12, height = 6)

